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Abstract 

We study energy localization on the oscillator-chain proposed by Peyrard and Bishop to model 
DNA. We search numerically for conditions with initial energy in a small subgroup of consecutive 
oscillators of a finite chain and such that the oscillation amplitude is small outside this subgroup 
on a long timescale. We use a localization criterion based on the information entropy and we verify 
numerically that such localized excitations exist when the nonlinear dynamics of the subgroup 
oscillates with a frequency inside the reactive band of the linear chain. We predict a mimium value 
for the Morse parameter (// > 2.25) (the only parameter of our normalized model), in agreement 
with the numerical calculations (an estimate for the biological value is fx = 6.3). For supercritical 
masses, we use canonical perturbation theory to expand the frequencies of the subgroup and we 
calculate an energy threshold in agreement with the numerical calculations. 
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I. INTRODUCTION 



A plethora of chemical processes involving the DNA macromolecule are known [1-3], 
for example the existence of denaturation bubbles containing a few broken H-bonds, and 
the transcription process triggered by the bonding of the biochemical complex to a specific 
region of the DNA (the so called TATA-box). The oscillator-chain model for DNA [4] 
was first proposed to study the thermal denaturation of the DNA macromolecule, i.e. the 
separation of the two strands. The dynamics of this model was first approximated with 
soliton techniques [5-8]. Our motivation for the present work was to study this model with 
methods of finite-dimensional dynamical systems, which could later be extended to a realistic 
model of the DNA without translation symmetry. In this work we consider a finite chain of 
N oscillators with initial condition restricted to a small group of n « N of consecutive 
oscillators. We define a localized motion as one in which the amplitude of oscillation is small 
outside a group of n max oscillators for all times, with n < n max « N , and we introduce 
a numerical criterion to quantify localization based on the information entropy. We use 
the correspondence conjecture (CC) of Flach et al. [9, 10] that the nonlinear dynamics 
of the isolated group of n oscillators must have frequencies inside the reactive band of 
the linearized chain for localization to be possible. Within this conjecture, we show that 
there is a minimum value for the Morse parameter (the only parameter of the model) for a 
localized excitation to be possible. The predicted value \i = 2.25 agrees with our numerical 
calculations. Last, since the linear frequencies of the isolated n-system lie in the dispersive 
band, an immediate consequence of the (CC) is that there must be a critical nonzero energy 
for localization (namely, for at least one of the n frequencies to exit the dispersive band). 
For supercritical values of the Morse parameter, we use canonical perturbation theory to 
evaluate the frequency shifts and we predict a threshold energy for localization in agreement 
with the numerical calculations. 

For the normalizations that follow, the most convenient way to introduce the Peyrard- 
Bishop (BP) model [4] for a DNA macromolecule is by the Lagrangian 

Lpb = 2j"2"("i + ** ? ) ~ 2^ Ui+1 ~ Ui ^ ~ 2^ m ~ Vi ^ ~ D i ex P a ( v i- u i) - !] 2 > (!) 

i=l 

where Ui and Vj denote the relative displacements of the nucleotidic bases at sites i and 
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j of each DNA strand, respectively. The number N denotes the number of sites in each 
strand of the DNA and can be as large as N ~ 10 9 . For practical reasons we perform 
our numerical experiments with up to N = 500. The masses of the bases have a common 
value m, the constant k is the longitudinal elastic constant and the parameters D and a 
define the Morse potential describing the transversal H-bonds linking the two chains. The 
experimental values for these parameters have been discussed in the literature[3, 13-15]: the 
mass of the base pairs is about (300 a.m.u. = 5.010~ 25 A;g), and the linear spring constant is 
(O.OAeV/ A 2 ). The hydrogen bond is modelled by the Morse potential with a = AA5A' 1 and 
for D we take an average of the value for the guanine-citosine basis (G-C) and the value for 
the thymine-adenine, A-T base pair, D = 0.04eV [13] 

By means of a rotation of coordinates defined by 

Xi = (ui + Vi)/y/2, (2) 
yi = {ui - Vi)/V2, 

the PB Lagrangian (1) splits into the sum L BP = L x + L Y , with 

N 

L X = Y, ™x 2 /2 - k{x l+1 - Xi ) 2 /2, (3) 
depending only on the x coordinates and with Ly depending only on the y coordinates as 

N 1 k 

Ly = J2 2 m ^ ~ 2 {Vl+l ~ Vi)2 ~ jD [ ex P(" v/ 2^) " X ] 2 ' ( 4 ) 
i=i 

Lagrangian Ly can be normalized by introducing a dimensionless time parameter r = 
a/ k/mt and a dimensionless coordinates ^ = V2ayi. The above scalings bring L Y of Eq. 
(4) to the normal form Ly = with 

^11 2 

i=i 

where dot denotes derivative with respect to r. Our normalization differs from that of [3], 
and it was chosen such that the quartic approximation to Eq. (4) has the form of the Klein- 
Gordon oscillator chain studied in [11, 12]. We henceforth study a chain of iV sites with 
periodic boundary conditions described by Lagrangian (5), a dynamical system depending 
on the single parameter /i 2 = ADa 2 /k (henceforth called the Morse parameter). Using the 

4 



values of the literature [13-15], we estimate a realistic biological value for /x to be /x = 6.3, 
and the scaling factors for the units of time and energy to be 0.810~ 12 s and 2.010~ 3 el / 
respectively. 

For small amplitudes, the normal mode frequency spectrum of (5) is [12] 

u(k) = \Jy? + 4sin 2 (fc/2), (6) 

where k = (jn/N), j = 1, ...N. The range of normal mode frequencies /x < uo < a/z^ 2 + 4 
constitutes the dispersive band, while the two relations < u < /x and > a//x 2 + 4 define 
the lower and upper reactive bands, respectively. One expects that localized motions of the 
chain with frequency components inside the dispersive band will give rise to quasi-normal 
modes excitations, which are typically delocalized in space. In such a way the localized state 
looses its energy in the form of radiation and spreads out. In contrast, localized excitations 
displaying only frequency components inside the two reactive bands are expected to preserve 
localization for long times. 

There is a large body of studies of one-dimensional chains, investigating the energy in- 
terchange among the "linearized " system modes [16-23]. For initial energy in a few low 
frequency modes, one of us (J. DL) has developed theoretical descriptions for energy spread- 
ing among modes, valid in various energy ranges, which were compared to numerical results 
for the FPU chain [16-18] and for the Klein-Gordon chain [12]. If the energy is initially 
placed in high frequency modes, the dynamics is transiently mediated by the formation of 
unstable nonlinear structures [20-23]. The mode energy is found to distribute itself first into 
a number of structures, localized in space, each consisting of a few oscillators, which coalesce 
over time into a single localized structure, a chaotic breather (CB). Over longer times the 
CB is found to break up, with energy transferred to lower frequency modes. Recently, there 
have been another set of studies of the discretized Klein-Gordon equation, from the perspec- 
tive of studying the stability of breathers, which are chosen as initial conditions [24, 25]. 
For a more comprehensive discussion of the extensive research on the dynamics of oscillator 
chains, see for example Ref. [12]. 

The study of soliton solutions of the nonlinear PDE's obtained by multiple-scale expan- 
sions constitutes at present the main line of study of the nonlinear dynamics of the DNA 
models[5, 6, 8]. Even though the use of modulation equations and soliton theory does 
furnish interesting results, we made here the choice to follow a different approach, based on 
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normal form methods for low-dimensional Hamiltonian dynamical systems. The reason for 
this choice is that, as is well known, multiple-scale expansions are only valid for initial data 
varying slowly in real space and quasi-monochromatic in Fourier space, while here we are 
interested in the evolution of initial excitations strongly localized in real space (delocalized in 
Fourier space). Moreover, the method used here displays the further advantage that it could 
be suitably extended to inhomogeneous chain models describing a realistic DNA molecule. 



II. THE CORRESPONDENCE CONJECTURE REVISITED 



In what follows we revisit the correspondence conjecture (CC) of Refs. [9, 10] in the 
light of canonical perturbation theory. With reference to system (5), let us initially displace 
from the equilibrium position £ = only a finite number n « N of consecutive particles. 
For such initial datum, instead of studying the dynamics of the full chain, involving a large 
number of degrees of freedom, we study the dynamics of the subsystem defined by the 
Lagrangian (5), with the sum restricted to the sites corresponding to the degrees of freedom 
initially excited, and with fixed end boundary conditions for the next neighbors. Such an n- 
degree of freedom subsystem is thought of as isolated and having energy E. This subsystem 
can be regarded as a perturbation of n linearly coupled oscillators, whose normal mode 
frequencies u>i, . . . , uj n are shown to lie inside the dispersive band of the larger iV-chain. For 
sufficiently low energies, the dynamics is quasi-linear and its frequency spectrum is close 
to the normal mode frequencies tvi, . . . ,u n . According to the (CC), if one uses such initial 
conditions for the larger lattice, the normal modes of the larger chain are excited and the 
initial excitation will spread out, which is what we observe numerically. With increasing 
energy, the effect of nonlinearity becomes prominent, and the frequency spectrum is modified. 
In the absence of resonances of third and fourth order (at least) in the harmonic spectrum 
of the subsystem, the modes preserve their identity and we can follow their frequency shifts 
inside the dispersive band. According to the (CC), one has localization for initial data at a 
given energy when the frequencies of the corresponding motion of the subsystem is outside 
the dispersive band. Of course this can happen only if the energy (i.e. the nonlinearity) is 
high enough. 

The conditions required on the frequency spectrum for multi-periodic oscillations are 
much more restrictive and the localization properties of such states can be very weak[10]. 
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For this reason, we restrict the analysis to periodic oscillations only, which amounts to 
look for periodic orbits of the subsystem whose frequency and harmonics lay outside the 
dispersive band. This analysis is detailed below. 



A. Analysis of the finite subsystems 



For the theoretical analysis we consider localized excitations where the amplitude of 
oscillation is small for sites outside a group of 1 + 2r modes, say \i\ > r (our i runs to both 
negative and positive directions and the central particle is % — 0). Under this conjecture the 
dynamics for sites on the right-hand side of the group {% > r) can be approximated by a 
linear chain driven by the given oscillation of oscillator r (while the same can be said of the 
oscillators on the left-hand side, i < —r). The equation of motion for the linearized chain 
can be derived from the Lagrangian (5) by expanding the exponential 

& = &+i + &-i- (2 + A* 2 )& , *>r (7) 

where the above linearization holds only for the oscillators outside the subgroup, which are 
supposed to oscillate with a small amplitude (i > r). The coordinate £ r (i) of oscillator r 
entering into Eq. (7) must be given a priori as a known forcing term. To solve Eq.(7) by 
Fourier transform we define the two-component vector 

X i+1 = ■ (8) 

It can be shown with the help of Eq. (7) that Xi+i sa tisfies the linear matrix iteration law 




Xi+i = \Xi, (9) 



where u 2 = 2 + fi 2 . For example in the case of a monochromatic forcing, Xr+ii 00 ) i s nonzero 
only at a single frequency u), and for the iteration of Eq.(9) to produce a bounded amplitude 
for sites of a large i it is necessary that (uo 2 — ou 2 ) 2 > 4, which is the definition of the reactive 
band (as opposed to the radiation band defined by [i < uj < a//i 2 + 4). If the forcing has 
several large Fourier components, the first large component might be in the lower reactive 
band {uj < /i ), while the other important harmonics could be in the upper reactive band 



U = TT 
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The first subsystem we consider here (henceforth called the 1-system) is defined by r = 
and consists of the nonlinear oscillation of a single particle of coordinate x a (t) with fixed 
ends = £ : = ). This nonlinear dynamics can be derived from Lagrangian (5) and it 
is also described by the Hamiltonian 

# =^ 2 + ^ + ^(exp(-U-l) 2 - (10) 

The frequency of oscillation for the periodic motion of Hamiltonian (10) can be determined 
by a simple quadrature for any energy, by the formula 

-, -l 

& ■ (ii) 

,„, ; „\/2B-2{;-^ (eX p ( _ { j_ 1) 2_ 

In figure 1 we plot this frequency as a function of the energy for several values of the 
parameter \x to illustrate that it is always inside the radiation band for fi < 2.25 at any 
energy. This is then the minimum value for the \i parameter where localization is possible, 
as predicted by the correspondence conjecture for the simple 1-system. It turns out that 
the biological value is /i — 6.3 > 2.25 , in agreement with this theory. Another agreement 
with this simple theory is discussed in the numerical section, as the numerical searches never 
found a localized state with \x < 2.5. 

For supercritical values of \x (// > 2.25 ), the frequency (11) is in the lower reactive 
band for a sufficient large energy. The frequency of small oscillations (zero energy) is easily 
obtained by expanding Hamiltonian (10) to quadratic order, and is u = v/2 + /i 2 > pi. The 
next correction for small energies can be obtained by expanding Hamiltonian (10) to fourth 
order in £ G as 



Hl 4) = \v 2 + \<j£ - ^ + 7 -£& (12) 

where the superscript and subscript on H refer to the order of the expansion and to the 
1-system, respectively. Introducing action-angle variables and using standard canonical 
perturbation theory [26], we find that the normal form of Hamiltonian (12) up to second 
order in the action variable J is 

H[ 4) =u J-g(fi)J 2 , (13) 

with 

- WTW- (14) 



Notice that for supercritical values of // (/i > 2.25), the coefficient g(/j,) as defined by Eq. 
(14) is positive (g(/i) > 0), such that the nonlinear frequency decreases with increasing 
energy Defining the nonlinear frequency by f2 = dH^ /dJ, the conditions = E c and 
f2 = dH^ /dJ = /i determine the minimum energy E c to be 

The interpretation of Eq. (15) is as follows: If the isolated nonlinear oscillator of the 1-system 
has an energy E > E c , its frequency is in the reactive band (Q < fi) and we expect that the 
corresponding type (i) initial condition should produce a localized excitation, according to 
(CC). This determination of the critical energy is compared to the numerical results in the 
following section, and it turns out to be short by a factor of two. The explanation for this is 
that the subsystem consisting of a single oscillator looses a significant amount of energy to 
the immediate neighbors, such that one could expect a higher critical energy. It turns out 
that the value E c ~ 1 predicted by Eq.(15) is precisely a factor of two short of the numerical 
value E c ~ 2 for any value of fi. Our simple theory is then seen to be only in qualitative 
agreement with the numerical calculations. A better approximation should be given by a 
subsystem consisting of three particles with fixed ends, which is our next subsystem. 

We consider another subsystem (henceforth called the 3-system), consisting of three os- 
cillators along the symmetric motion defined by £_ x = £ v The Lagrangian equations of 
motion derived from (5) with the condition = £ x correspond to the following two-degree 
of freedom Hamiltonian 

Hs = Ipl + lpl + Zl + iL-Zi) 2 (16) 
V(exp(-ei) - If + y(exp(-0 - l) 2 . 
The two frequencies of the quasi-periodic linear motion at zero energy are 

u 1 = ^^ + 2-V2, (17) 
tu 2 = \Jy? + 2 + V2, 

which are inside the dispersive band for any \i. 
For small energies subsystem (16) is a perturbation of two harmonic oscillators with 
frequencies c^iand u 2 inside the dispersive band of the whole linearized chain. To compute 
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the leading contribution to the frequency shift of each oscillator we must evaluate the next 
frequency correction in powers of the actions. One can check that there is no resonance up 
to fourth order involving the linear part of (16), i.e. LO\juj 2 ^ 1/2,1/3. Using canonical 
perturbation theory [26] we can remove the cubic term from (16), average the quartic term 
and express the normal form of Hamiltonian (16) up to second order in the actions as 

#3 4) (Jl, J 2 ) = LO\J\ + UJ 2 J 2 ~ Cl J\ ~ C 2 Jl ~ C12 Jl J2, (18) 

where the J's are the action variables and ci, c 2 and C12 are given by 



C12 



c 2 = 



3fi 2 [12/i 6 + (65 - 26-s/2)/ + (73 - 75y/2)(i 2 - (42 - 35^)] 

3/i 2 [36// + 147/i 6 + 4/i 4 - 278// 2 - 98] 
16ujIuj 3 2 {Auj 2 - J*) (4^1 - ujj) ' 
3fi 2 [12/i 6 + (65 + 26V2> 4 + (73 + 75y/2)f? - (42 + 35^2)] 



64ujul (4cj 2 - uj) 
The nonlinearly modified frequencies are given by 



dH 

Q 1 = — — = LO X - ICxJx - C 12 J2, (20) 
OJl 

dH 

^2 = TTT = ^2 - 2C 2 J 2 - C12J1. 

For supercritical values of fi the coefficients of (19) are all positive, such that the fre- 
quencies of (20) are decreasing functions of the energy. The two periodic orbits branching 
from the linear modes of (18) are obtained by setting one of the actions of (18) to zero. For 
example by substituting J 2 = into (18) we obtain 

Hi 4) (J 1 ,0)=uj 1 J 1 -c 1 Jl (21) 

such that the critical energy predicted by setting = /i is 

E * = adfJ* = iz:/!. (22) 

4ci 4ci 

For the other periodic orbit we substitute Ji = into (18), yielding 

i^ 4) (0, J 2 ) = uj 2 J 2 -c 2 J 2 , (23) 
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and the critical energy predicted by setting f2 2 = [i is 

E ( 2 ) = = 2W2. (24) 

4c 2 4c 2 

For values of /i in the interval (2.5 < fi < 30) one sees that the values of c\ and c 2 are close 
to the limiting values c\ ~ c 2 ~ ^, while ci 2 has the limiting value ci 2 ~ |. (It is easy to 
obtain this limit by setting uj\ ~ uo 2 ~ /i and 4a;^ — cj^ ~ 4u;2 — c^i ~ 3/i 2 into the formulas 
of (19) ). The limiting values for the critical energies are E { c ] = 4(2 - y/2)/3 ~ 0.78 and 
E ( c ] = 4(2 + v^)/3 = 4.55. The critical energy E^ = 0.78, obtained for localized excitations 
generated by J 2 = 0, agrees within twenty- five percent with the numerical calculations of the 
next section, which determine E c ~ 0.6. For initial conditions in the 3-system, the energy 
leaking out is compensated by a negative interaction energy of the 3-system with the rest, 
such that the energy inside the 3-system is actually larger than the total energy (this explains 
how we have overestimated the critical energy). The reason for this better agreement is still 
that, by increasing the subsystem size, the interaction energy with the immediate neighbors 
(whatever its sign) becomes less important. Models with more oscillators in the subgroup 
should furnish even better approximations, but they are harder to work out analytically 
and the corresponding type(n) initial conditions are computationally more expensive to 
investigate. 

III. NUMERICAL RESULTS 

We present numerical results for the DNA oscillator chain, with initial condition in two 
different types of oscillator groups. All of our numerical integrations were performed with a 
tenth order symplectic Runge-Kutta-Nystrom integrator [27]. The high-order integrator can 
take very large steps, of about 0.6 of the shortest linear period and still conserves energy 
with a precision of 10~ 10 even after integration times of 10 10 . 

A. Macroscopic quantities 

The dynamics of the full chain described by Lagrangian (5) is described by the following 
Hamiltonian 

N/2 1 1 1 
H = E 2^ + 2 fe+1 _ Q2 + ^^'^ " lf ■ (25) 

i=-N/ 

11 



In numerical experiments the instantaneous values of the on-site oscillator energies Ei, 
i = 1, N, is usually calculated as 

Ei = \p1 + - + l^i+i " + l^eM-Q - I) 2 , (26) 

where we include fifty percent of the interaction with the oscillator at each side, such that 
the sum of the Ei is the constant total energy. Over short times the instantaneous and 
average values are nearly the same. The information entropy is defined by 

N 

S = -J2e l lne l , (27) 
i=i 

where = Ei are the normalized instantaneous oscillator energies. In a typical 

situation where the total energy is distributed among r < N oscillators, r of the are of 
order 1/r and the remaining are negligible, such that Eq. (27) predicts S ~ ln(r). This 
motivates the definition of 

N osc = exp(S), (28) 

as the effective number of oscillators sharing the energy. It is also convenient to define the 
normalized parameter 

n sc = N osc /N. (29) 

The normalized parameter n osc varies from to 1, because the entropy of Eq.(27) is al- 
ways less than ln(iV). The instantaneous value of n osc does not asymptote to one, due 
to fluctuations. To calculate the effect of fluctuations we introduce a deviation <5ej from 
equipartition = e + bt{. Substituting this into (29), expanding the logarithm function as 
ln(l + Sei/e) = Sei/e — (l/2)(5ej/e) 2 and performing the summation over i yields 

n osc = ^exp{-iVeln(e) - iV(5e) 2 /(2e)} = exp{-iV(5e) 2 /(2e)}. (30) 

Taking e — l/N and making the assumption of normal statistics, that for each normal mode 
(5e) 2 = e 2 (which is true only for the linearized lattice dynamics ), we see that iV cancels 
giving an asymptotic value n osc = exp(— 0.5) = 0.61. This calculation shows that the result 
does not depend on the number of oscillators if N is large and also shows why the value is 
different from unity. More accurate calculations have been made, including the nonlinear 
terms in the oscillator calculation, yielding[23] 

n osc = 0.74 (31) 

12 



at equipartition of energy among the oscillators. These values have been checked numerically, 
giving good agreement [23]. 

Numerical experiments show that for a randomly chosen localized initial condition, the 
value of n osc usually starts to increase and reaches the equipartition value n osc = 0.61 in a 
time of the order of iV , which is the typical spread time. Our localization criterion is that 
a state is localized when the value of n osc is significantly less than the equipartition value 
0.61 for more than 50iV periods Tf of the fastest linear mode (T/ = 2tt/ a//i 2 + 4 ~ 1). For 
the computationally accessible finite values of N (of the order of 100), the smallest value 
of n osc is obtained for localization in a single site, n osc = 1/N, which is an extreme value. 
Given that n osc = 0.61 means equipartition, our practical criterion is n osc < rj L = (20/N) for 

t < 50N. With this criterion we give the state some room to breathe, allowing the energy 
to spread over twenty oscillators and then to shrink again to a smaller average value. This 
practical criterion excludes either states localized in more than twenty oscillators or states 
that would have sudden derealization bursts, which was never observed numerically. In 
the numerical calculations we use a logarithmic scale for the increasing time, in the natural 
units of Lagrangian (5). The rapid fluctuations of the instantaneous values are smoothed 
by taking the average of the last five instantaneous values of n osc , which are evaluated at a 
rate of 25 points per decade in time (at every integer value of 251n(t)). 

Our numerical experiments integrate the dynamics of Lagrangian (5) for a chain of N 
oscillators with periodic boundary conditions and we shall use two types of initial condi- 
tions, defined as follows : (i) Initial conditions produced by giving a nonzero position and 
momentum to a single oscillator and a null value for the positions and momenta of all other 
oscillators (the value of N osc at t — is 1). (ii) Symmetric initial conditions produced by 
giving a nonzero value for three consecutive oscillators with the symmetry x-i = x\ and 
p_i = pi (the value of N osc at t — is 3). For example we have used r] L = 0.2 and we 
started several (about 50) initial conditions of type (i) with a given energy. For each initial 
condition we calculate n osc along the numerical integration and we stop the integration at 
the first time that n osc becomes larger than rj L = 0.2, defining a derealization time for that 
initial condition. The maximum value of the derealization time ( T max ) among the 50 
initial conditions of the same energy is our measure of localization. Typically, for a chain of 
N = 100 oscillators, this value is about T max ~ 100 = N for subcritical energies, than there 
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is a rapid transition where this value climbs to above T max = 5000. In practice, it is nec- 
essary to stop the numerical integration in the supercritical region whenever T max reaches 
a maximum value, and we have used T max = 53iV as a good computationally accessible 
large number (53iV = 5300 if N — 100). We experienced with a much higher threshold for 
T max of about lOOOiV and obtained virtually the same type of transition, but the numerical 
experiment becomes very time consuming. The question if this localization time is either 
infinite, exponentially long, or simply a very large value is not addressed in the present 
work. We have also varied the threshold value of n osc among the values r/ L = 0.15, r\ L = 0.2 
and rj L = 0.25 and we obtained the same transition line. We used for N the three values 
N = 100, N = 200 and N = 500 and obtained virtually the same transition lines for /i > 3. 
A comprehensive statistical analysis has not been performed due to the very long times for 
some runs. Spot checks for a few cases indicate that the spread from varying iV and i] L is 
less than some few percent if // > 3 and r\ L < (25/iV). For the region 2.5 < /i < 3 there 
can be significant changes in the critical energies determined by the above procedure. This 
is because close to the critical value \i = 2.25 the localization length becomes long, and 
in a lattice with a small iV this localization is confused with equipartition by our criterion. 
Interesting to recall that biology chose the safe value of \x = 6.3 possibly for the same reason. 

In Fig. 2 we plot the value of T max < 53iV among 49 initial conditions of type (i) as a 
function of the energy for \i = 6.3. Notice the pronounced jump in T max which is a signature 
of localization. We define the critical energy by the inflection point of the T max curve, which 
from Fig. 2 is about E ~ 2.3. This same discontinuous behavior of T max is observed in the 
numerical calculations for 2.5 < /x < 30, and in Fig. 3 we plot the critical energy determined 
by the inflection point of T max and the theoretical predictions for the 1-system, Eq.(15), 
versus \x . The numerically determined critical energy is about twice the predicted for the 
simple 1-system by perturbation theory. This effect is due to the fact that for type (i) initial 
conditions a substantial part of the energy leaks to the immediate neighbors even when there 
is localization, such that the total energy of the system at localization is significantly larger 
than the energy of the 1-system. 

In Fig. 4 we plot T max < 53iV among 49 type (ii) initial conditions, as a function of 
the energy for \i = 3.0, and N = 100, illustrating the same jump that is our signature of 
localization. The critical energy predicted by the inflection point of Fig. 4 is E = 0.75. 
In Fig. 5 we plot the critical energy determined by the inflection point of T max and the 
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theoretical predictions of the 3-system versus /i. The theoretical prediction for the 3-system 
agrees with the numerical results within twenty five percent. The approximation is better 
than in the case of the 1-system because less energy leaks out of the 3-system. For type(ii) 
initial conditions there is an interaction term in the total energy that increases the energy 
inside the 3-system above the total energy, but the predicted energy is now only twenty five 
percent wrong. 

In Fig. 6 we plot the modulus of the complex Fourier transform of the coordinate of the 
central oscillator for an initial condition of type (i) of a lattice with /j, — 6.3 , N — 100 and a 
subcritical energy E = 0.1. Notice that the Fourier transform develops nonzero components 
inside the conduction band 6.3 < uo < 6.61, as illustrated in the insert to Fig. 6. 

In Fig. 7 we plot the Fourier transform of the coordinate of the central oscillator for a 
localized initial condition of type (i) in a lattice with /i = 6.3, N = 100 and a supercritical 
energy E = 3.0, which has a primary peak at u = 6.0 < /i and goes to zero already at 
uj = 6.25 < /i, in accordance with the (CC) conjecture. 

Last, in Fig. 8 we plot the surface of section of the 3-system with n = 6.3 at the 
supercritical energy E = 3, showing very little stochasticity, to illustrate that localization 
has nothing to do with stochasticity within the subgroup, as discussed in Ref. [10]. 

IV. DISCUSSIONS AND CONCLUSION 

At a supercritical energy, by searching among 49 initial conditions of type (i), for example, 
we have found several initial conditions that stay localized for more than 10 5 natural units. 
Using a numerical search that varies the initial condition in the neighborhood of an original 
localized condition [28], in a way that maximizes the localization time, we could easily find 
other initial conditions that stays localized for a much larger time, of the order of 10 7 . These 
refined initial conditions become restricted to narrow domains, and we believe that the study 
of timescales for a localized excitation in a chain with a finite N should start from here in 
future work, for example to test if one can increase this time arbitrarily. 

At a finite N, if the energy of a type (ii) initial condition does not localize in the original 
3-system, it will leak out to the other (N/3) 3-systems of the chain. A simple condition 
for these other 3-systems to be "sufficiently linear" is then that the total energy be less 
than NE c /3 (such that the other 3-systems display a quasi-linear motion). This intensive 
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condition, E < NE c /3, is important to remember in numerical experiments with a finite 
lattice. For example for a chain of N = 100 oscillators, this means E < 26.4 and in our 
numerical experiments we have always stayed well below this energy. 

The critical value of fi for localization (// = 2.25 ) is in agreement with the numerical 
calculations, as we never found localization below \i = 2.5. In the region 2.5 < \i < 3.0, the 
numerical results indicate that the localization length is very large, which requires numerical 
experiments with large values of N. The value \x = 6.3 estimated from the biological 
measurements is far from the critical and in a region where localization length is small, such 
that we predict a robust localization from the above DNA model. The threshold energy 
for localization at fi = 6.3 is 2.2 units or 4.4 x 10~ 3 eV ( 0.17/^T at room temperature). 
This means that localization is possible at room temperature, as predicted by our model. 
Last, the localization time found numerically is greater than 10 5 time units or 10~ 7 seconds, 
enough for the biochemical mechanisms to operate. Such localization can be related to the 
bubbles in DNA and it would be an auxiliary mechanism in the transcription process. 
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VI. FIGURE CAPTIONS 

Fig. 1 Frequency of the nonlinear 1-system divided by fi, (co/fj), plotted as a function of 
the energy for fi = 2 ( dashed line), fi = 2.25 (solid line) and fi = 2.5 (dotted line ). The 
horizontal solid line is cu/fi — 1. Notice that at fi — 2.25 the frequency line is only tangent 
to the critical line. Arbitrary units. 

Fig. 2 T max as a function of the energy for type (i) initial conditions at N = 100 and 
fi = 6.3. Arbitrary units. The squares represent numerical calculations and the solid line is 
spline interpolation. 
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Fig. 3 Numerically calculated critical energy for type (i) initial conditions with N = 100 
(triangles), N = 200 (stars) and N = 500 (circles) as a function of the energy. Also plotted 
is the critical energy predicted by the 1-system (squares ). Arbitrary units. 

Fig 4. T max as a function of the energy for type (ii) initial conditions at N = 100 and 
\i = 3.0(triangles). Arbitrary units. The triangles represent the numerical calculations and 
the solid line is spline interpolation. 

Fig. 5 Numerically calculated critical energy for type (ii) initial conditions with N = 100 
as a function of fi (triangles) and critical energy predicted by the 3-system E^ (circles) , as 
a function of \i. Arbitrary units. 

Fig. 6 Fourier transform of an initial condition of type (i) with E — 0.1 (subcritical) for 
a lattice with \x = 6.3 and N = 100. Plotted is the modulus F(u) of the complex Fourier 
transform. The insert magnifies the region near uo = 6.3 to display that F(uo) is not zero 
inside the conduction band. 

Fig.7 Fourier transform of an initial condition of type (i) with E = 3.0 for a lattice with 
pi = 6.3 and N = 100. Plotted is the modulus F(uo) of the complex Fourier transform. The 
insert magnifies the region near u = 6.0 to display the peak of F(uo) at uo = 6.0 < /i. Notice 
that F{u) vanishes above uo = 6.25 < //. 

Fig. 8 Surface of section of the symmetric 3-system at /i — 6.3 and E = 3.0, showing 
little stochasticity. Arbitrary units. 
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